setwd("/home/jar/zxjar/paper_helper/userdata/0/78")
library(readr)

library(haven)
library(plyr)
library(dplyr)
#加权分析用
library(survey)
#tableby
library(arsenal)
library(gtsummary)

#setwd('F:/Rproject/Adultnon-alcoholic')

paper.data<- read_csv("fromdata.csv")

NHANES_design <- svydesign(
data = paper.data,
ids = ~SDMVPSU,
strata = ~SDMVSTRA,
nest = TRUE,
weights = ~WTINT2YR,survey.lonely.psu = "adjust"
)



tbl_svysummary(NHANES_design,
              missing = 'no', 
              # 通过属性分组
              by = fibrosis_group,
              # 分区yes 需要YES 和NO 如果采用Yes 和No 不知道是不是有变量冲突导致分类不显示 只显示全部的
              # 进行设置显示维度
              include = c(RIDAGEYR,Sex,Raceethnicity,Educationlevel,BMIgroup,rpa,smoke_group,diabetesindex,
                          LBXVIC,LBXSATSI,LBDHDD,LUXCAPM),
              # 批量给变量起别名
              label = list(RIDAGEYR ~ "Age (years)",
                           LUXCAPM ~ "Median CAP (dB/m)",
                           Sex ~ "Gender, n (%)a",
                           rpa ~ "Recreational physical activity, n (%)a",
                           Raceethnicity ~ "Race/ethnicity, n (%)a",
                           Educationlevel ~ "Education level, n (%)a",
                           BMIgroup ~ "BMI group, n (%)a",
                           smoke_group ~ "Smoking status, n (%)a",
                           LBXVIC ~ "Serum vitamin C (mg/dL)",
                           diabetesindex ~ "Diabetes, n (%)a",
                           LBXSATSI ~ "ALT (IU/L)",
                           LBDHDD ~ "HDL-cholesterol (mg/dL)"),
              statistic = list(
                # 分别对应数值型和分类变量 后是需要显示表达式
                all_continuous() ~ "{mean}±{sd}",
                # 样本是加权倒置n也是加权的 下面参数是显示
                all_categorical() ~ "{n_unweighted} ({p}%)"
              ),
              digits = list(all_continuous() ~ 2
                            # ,all_categorical() ~ 3
                            )
              )%>%
   add_n(
     statistic = "{N_nonmiss_unweighted}",
     col_label = "**N**"
   
   ) %>% # 添加非NA观测值个数
  # add_p() %>%
  # 这里p值不能用p 根据原文说的是加权得来的 需要用 tbl_svysummay进行 这里除了p值 其他是正确的
  add_p() %>% # 添加P值
  add_overall() %>%
  #modify_header(all_stat_cols() ~ "**level**, N = {n_unweighted} ({stype_percent(p)}%)" ) %>%
  # 由原来的加权显示改成非加权显示原始值 
  modify_header(stat_1 ~ "**No** (n={n_unweighted})") %>%
  modify_header(stat_2 ~ "**Yes** (n={n_unweighted})") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Significant fibrosis**")%>%as_flex_table()%>%flextable::save_as_html(path = '2.9tab1.html')
  # 编辑表头测试
  #modify_spanning_header(stat_1 ~ NA, update = all_stat_cols() ~ "**Significant fibrosis**") %>%
  #modify_footnote(update = all_stat_cols() ~ "测试编辑表格下面值")
  #%>%
  #bold_labels()  #label 粗体

# 线性回归模型进行测试
# Serum vitamin C
# NHANES_designfg1 <-subset(NHANES_design,fibrosis_groupint==1)
# model1
glmsvc <- svyglm (fibrosis_groupint ~ LBXVIC+LBXVIC.quantile.var ,design = NHANES_design,family=quasibinomial )
model1 <-tbl_regression(glmsvc,exponentiate = T,label = list(LBXVIC ~ 'Serum vitamin C (mg/dL)',LBXVIC.quantile.var ~ 'Serum vitamin C (mg/dL, quartile)') )%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# model2

#  模型2 年龄性别和种族
# RIDAGEYR +  Sex + Raceethnicity

# 多元由原来的分类变量 试试连续型变量
glmsvcmodel2 <- svyglm (fibrosis_groupint ~ LBXVIC+LBXVIC.quantile.var+RIDAGEYR +  Sex + Raceethnicity ,design =NHANES_design,family=quasibinomial )

model2<-tbl_regression(glmsvcmodel2,include = c(LBXVIC,LBXVIC.quantile.var),label = list(LBXVIC ~ 'Serum vitamin C (mg/dL)',LBXVIC.quantile.var ~ 'Serum vitamin C (mg/dL, quartile)') ,exponentiate = T)%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


# model3 
#模型3  年龄、性别、种族/族裔、BMI、教育水平、吸烟状况、休闲体育活动、高密度脂蛋白胆固醇和糖尿病
# RIDAGEYR +  Sex + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index

glmsvcmodel3 <- svyglm(fibrosis_groupint ~ LBXVIC+LBXVIC.quantile.var+RIDAGEYR +
                          Sex + Raceethnicity + BMIgroup + Educationlevel +
                          smoke_group+ rpa+  LBDHDD + diabetes.index ,design = NHANES_design,family=quasibinomial)

model3<-tbl_regression(glmsvcmodel3,include = c(LBXVIC,LBXVIC.quantile.var),label = list(LBXVIC ~ 'Serum vitamin C (mg/dL)',LBXVIC.quantile.var ~ 'Serum vitamin C (mg/dL, quartile)') ,exponentiate = T)%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


alllei2<-tbl_merge(tbls = list(model1,model2,model3),tab_spanner = c('Model 1OR (95% CI)','Model 2 OR (95% CI)','Model 3OR (95% CI)'))
alllei2%>%as_flex_table()%>%flextable::save_as_html(path = "2.9table2.html")

#model1 
# model2
# RIDAGEYR +  Sex + Raceethnicity
# model3 
#模型3  年龄、性别、种族/族裔、BMI、教育水平、吸烟状况、休闲体育活动、高密度脂蛋白胆固醇和糖尿病
# RIDAGEYR +  Sex + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
#model1  这里面有个问题 如果你设置的截取数据在你分析协变量的维度里面是不能是已经拆分出固定的值 就比如数据里面性别选择了男 那么协变量不能有性别维度进行分析会报错

glmsvcmodel3male <- svyglm(fibrosis_groupint ~ LBXVIC
                           ,design = subset(NHANES_design,Sex =='male'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'male'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


glmsvcmodelfemale <- svyglm(fibrosis_groupint ~ LBXVIC,
                            design = subset(NHANES_design,Sex =='female')) %>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),
                 label = list(LBXVIC ~ 'female'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Overweight

glmsvcmodel3obese <- svyglm(fibrosis_groupint ~ LBXVIC
                            ,design = subset(NHANES_design,BMIgroup =='>=30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'obese'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3underweight <- svyglm(fibrosis_groupint ~ LBXVIC
                                  ,design = subset(NHANES_design,BMIgroup =='<25'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'Under/normal weight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3Overweight <- svyglm(fibrosis_groupint ~ LBXVIC
                                 ,design = subset(NHANES_design,BMIgroup =='25-30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'overweight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

model1<- tbl_stack(tbls = list(glmsvcmodel3male,glmsvcmodelfemale,glmsvcmodel3underweight,glmsvcmodel3Overweight,glmsvcmodel3obese),group_header = c('Stratified by gender','Stratified by gender','Stratified by BMI','Stratified by BMI','Stratified by BMI'))

#model2 
# RIDAGEYR +  Sex + Raceethnicity

glmsvcmodel3male <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR + Raceethnicity
                           ,design = subset(NHANES_design,Sex =='male'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'male'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


glmsvcmodelfemale <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR  + Raceethnicity,
                            design = subset(NHANES_design,Sex =='female')) %>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),
                 label = list(LBXVIC ~ 'female'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Overweight

glmsvcmodel3obese <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity
                            ,design = subset(NHANES_design,BMIgroup =='>=30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'obese'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3underweight <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity
                                  ,design = subset(NHANES_design,BMIgroup =='<25'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'Under/normal weight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3Overweight <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity
                                 ,design = subset(NHANES_design,BMIgroup =='25-30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'overweight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")



model2<- tbl_stack(tbls = list(glmsvcmodel3male,glmsvcmodelfemale,glmsvcmodel3underweight,glmsvcmodel3Overweight,glmsvcmodel3obese),group_header = c('Stratified by gender','Stratified by gender','Stratified by BMI','Stratified by BMI','Stratified by BMI'))


# model3 
# RIDAGEYR +  Sex + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
glmsvcmodel3male <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR  + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
                           ,design = subset(NHANES_design,Sex =='male'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'male'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


glmsvcmodelfemale <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR  + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index,
                            design = subset(NHANES_design,Sex =='female')) %>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),
                 label = list(LBXVIC ~ 'female'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Overweight

glmsvcmodel3obese <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity  + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
                            ,design = subset(NHANES_design,BMIgroup =='>=30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'obese'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3underweight <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity  + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
                                  ,design = subset(NHANES_design,BMIgroup =='<25'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'Under/normal weight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3Overweight <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity  + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
                                 ,design = subset(NHANES_design,BMIgroup =='25-30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'overweight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")



model3<- tbl_stack(tbls = list(glmsvcmodel3male,glmsvcmodelfemale,glmsvcmodel3underweight,glmsvcmodel3Overweight,glmsvcmodel3obese),group_header = c('Stratified by gender','Stratified by gender','Stratified by BMI','Stratified by BMI','Stratified by BMI'))



alllei3<-tbl_merge(tbls = list(model1,model2,model3),tab_spanner = c('Model 1OR (95% CI)','Model 2 OR (95% CI)','Model 3OR (95% CI)'))
alllei3%>%as_flex_table()%>%flextable::save_as_html(path = "2.9table3.html")






